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The evolution of a nonaxisymmetric bar-mode perturbation of rapidly rotating stars due to a secular 
instability induced by gravitational wave emission is studied in post-Newtonian simulations taking 
into account gravitational radiation reaction. A polytropic equation of state with the polytropic 
index n = 1 is adopted. The ratio of the rotational kinetic energy to the gravitational potential 
■ energy T/\W\ is chosen in the range between 0.2 and 0.26. Numerical simulations were performed 

' until the perturbation grows to the nonlinear regime, and illustrate that the outcome after the secular 

, instability sets in is an ellipsoidal star of a moderately large ellipticity ^ 0.7. A rapidly rotating 

(N . protoneutron star may form such an ellipsoid, which is a candidate for strong emitter of gravitational 

waves for ground-based laser interferometric detectors. A possibility that effects of magnetic fields 
^ , neglected in this work may modify the growth of the secular instability is also mentioned. 
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I. INTRODUCTION 



. Rapidly rotating stars are subject to nonaxisymmetric rotational instabilities. An exact treatment of these instabil- 
OO ' ities exists for incompressible and rigidly rotating equilibrium fluids in Newtonian gravity [1,2]. For these configura- 
, tions, global rotational instabilities arise from nonradial toroidal modes e™''' (m = 1, 2, . . .) when (3 = T /\W\ exceeds 
' a certain critical value. Here (p is the azimuthal coordinate and T and W are the rotational kinetic and gravitational 
potential energies (see Sec. II B for definition). There exist two different types of mechanisms and corresponding time 
^ ■ scales for the instabilities. One is the secular instability which, with a given value of m, sets in for a value of (3 larger 
O-i than a critical value (is and can grow in the presence of some dissipative mechanism, like viscosity or gravitational 
Q . radiation. The growth time is determined by the dissipative time scale, which is usually longer than the dynamical 
^ ' time scale of the system. By contrast, a dynamical instability, which sets in for a value of (3 larger than a critical 
c/3 value I3d(> Ps), is independent of any dissipative mechanisms, and the growth time is determined typically by the 

■ hydrodynamical time scale of the system. In this paper, we study the growth of a secularly unstable bar-mode for 
compressible stars whose growth is induced by the gravitational radiation. Hereafter, we only focus on the m = 2 
bar-mode since it is the fastest growing mode in most of rapidly rotating stars. (With regard to counter-examples in 

. ^ a narrow parameter range, see [3,4] about secular instability.) 

■ The criterion for the onset of the secular bar-mode instability for compressible rotating stars has been extensively 
determined by linear perturbative analyses in Newtonian theory [5-12], in post-Newtonian approximation [13], and in 
general relativity [14,15]. These studies indicate that the value of Ps for a bar-mode is « 0.14 in many cases, although 
it varies depending on the rotational law, equations of state, and effect of general relativity. In particular, the value 
of Ps could be decreased significantly below 0.14 for highly differentially rotating stars [12] and for a compact general 
relativistic star [14]. However, it is not increased much beyond 0.14 to our knowledge. Thus, in the presence of 
gravitational radiation reaction, the bar-mode perturbation can grow for rapidly rotating and secularly unstable stars 
with P > 0.14. 

The secular instability by gravitational wave emission may be relevant for rapidly and differentially rotating pro- 
toneutron stars formed soon after supernova collapse [4]. To study the growth of a secularly unstable bar- mode 
perturbation to a nonlinear perturbative regime, a longterm numerical simulation is necessary. However, such study 
has not been performed until quite recently even in the Newtonian or post-Newtonian gravity. (But, see [16] for a 
study of the incompressible case in which the basic equations reduce to ordinary differential equations, and [17] for a 
simulation of r mode instability.) 

There are several evolutionary paths which may lead to the formation of rapidly and differentially rotating neutron 
stars with /? > 0.14. Assuming mass and angular momentum conservation, /3 increases approximately as during 
stellar collapse where R denotes the radius of the collapsing core. In rotating supernova collapse, the core contracts 
from R ^ 2000 km to a few 10 km (e.g., [18-21] for Newtonian simulations and [23,24] for general relativistic 
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simulations), and hence, (3 may increase by two orders of magnitude. In reality, a large fraction of the mass is ejected 
in the supernova explosion, and hence, the relation between /3 and R will not be as simple as expected (e.g., [23]). 
Nevertheless, the value of /3 in the inner region that forms the protoneutron star should increase. The radius of 
a rotating progenitor star that forms a protoneutron star of radius ~ 10 km will be ^ 100 km, and hence, it is 
reasonable to expect that the value of (3 will increase at least by one order of magnitude. Thus, a stellar core with the 
effective value oi (3 'i> 10~^, which can be easily reached in differentially rotating cases [23], may yield rapidly rotating 
neutron stars which may reach the onset of secular or dynamical instability. Similar arguments hold for accretion 
induced collapse of white dwarfs to neutron stars [25] and for the merger of binary white dwarfs to neutron stars. 

Only if the value of /? is in a range between f3s and /J^, a rapidly rotating neutron star is secularly unstable but 
dynamically stable. Thus, it may be more likely that a rapidly rotating protoneutron star is dynamically unstable. 
The fate of dynamically unstable rotating stars against a bar-mode perturbation has been studied widely so far not 
only in Newtonian theory [26-30] but also in general relativity [31]. These studies have indicated that the value of (3d 
is ~ 0.27 in many cases except for highly differentially rotating stars [30] and that after the growth of the bar-mode 
perturbation, spiral arms are formed, and then, angular momentum is transported outward with mass-shedding. As a 
result, the angular momentum of around the center of the rotating star is decreased below the dynamically unstable 
limit to be a slightly nonaxisymmetric star in which the value of (3 is smaller than f3d- However, (3 remains still much 
larger than /3s. Therefore, the relic of a dynamical instability can still be subject to a secular instability. 

In this paper, we report the results of numerical simulations for the secular bar-mode instability of rapidly and 
differentially rotating protoneutron stars induced by gravitational radiation. The polytropic equation of state with 
the polytropic index n = 1 is adopted to approximately model a protoneutron star. The hydrodynamic simulations 
are performed in a (0-1-2.5) post-Newtonian framework; i.e., Newtonian gravity and gravitational radiation reaction 
are taken into account [32] . While protoneutron stars are fairly compact objects with reasonably strong gravitational 
field, Newtonian theory is expected to describe them up to errors of the order GAl/Rc? ^ 0.1-0.2 where G, c, M, 
and R are the gravitational constant, the speed of light, and mass and radius of the protoneutron stars, respectively. 

Although a longterm, stable, and accurate numerical simulations are now feasible (e.g., [35]), a simulation for the 
growth of a nonaxisymmetric secular instability in full general relativity is still formidable because of the following 
reasons, (i) it is required to take a large computational region that is extended to a wave zone to correctly compute 
gravitational radiation reaction which is the key process in this issue, and (ii) it is required to perform a very 
longterm simulation because the growth time scale of a secularly unstable perturbation is much longer (several orders 
of magnitude longer) than the dynamical time scale of the system. These two facts imply that computational resources 
much larger than the present best ones are necessary to stiidy this problem in fully general relativistic simulations. 
On the other hand, in a post-Newtonian simulation, (a) we can take into account the gravitational radiation reaction 
simply adding a radiation reaction potential, and hence, we do not have to take a large computational domain that is 
extended to a wave zone. Furthermore, (b) the magnitude of radiation reaction force can be artificially increased to 
shorten the radiation reaction time scale as short as a dynamical time scale. Specifically, we increase the magnitude 
of the radiation reaction term by a factor ^ 10 100 in this paper to accelerate the growth of a secularly unstable 
perturbation. (In the (0-1-2.5) post-Newtonian framework, increasing the radiation reaction force is equivalent to 
increasing the compactness of a star GM/Rc^; see Sec. II A.) The facts (a) and (b) enable to reduce the computational 
time significantly, and hence, the study of the secular instability for compressible stars is feasible in a post-Newtonian 
framework with current computational resources. 

The secular instability of rotating stars is also induced by viscous effects in a different way from that by gravitational 
radiation [1,36,4]: During the growth of the secular instability induced by gravitational radiation, vorticity is conserved 
but angular momentum is dissipated. On the other hand, in the case of viscosity, the angular momentum is conserved 
but the vorticity is dissipated. If the viscous dissipation time scale is as short as that of gravitational radiation, 
two effects interfere each other and suppress the growth of the secular instabilities [36]. An estimate suggests [4] 
that the viscous time scale in rapidly rotating protoneutron stars may not be as short as the dissipation time scale 
by gravitational waves, although it is not very clear. In this paper, we ignore the viscous effect as a first step and 
focus on the secular instability induced by gravitational radiation. Simulation adding the viscous term (solving the 
Navier-Stokes equations) is an interesting subject remained in the future. 

The paper is organized as follows. In Sec;. II, we describe basic equations, initial conditions, and the methods in 
numerical analysis. In Sec. Ill, the numerical results are presented paying attention to the growth time of a secularly 
unstable perturbation, the outcome after the growth of the secular instability, and gravitational waves. Section IV is 
devoted to a summary. Throughout this paper, we use the geometrical units of G = c = 1. 
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II. METHOD 



A. Basic equations 



Numerical simulations are performed in a (0+2.5) post-Newtonian framework [32]. Here, "0" implies "Newtonian" 
order and "2.5" the lowest order radiation reaction effect due to mass-quadrupolc gravitational wave emission. We 
adopt the following basic equations in Cartesian coordinates [32-34]: 

^ + ^=0, (1) 

dpuj d{puiv^ +PSi) _ _ djip + tpn) 

dt + dxj ~ ^ dxi ' ^ ' 

dpe dipe + Py ^dP djiP + V>fl) 

-m + dx^ = . (3) 

The first, second, and third equations are the continuity, Euler, and energy equations, respectively. Here, p, P, e, u', 
and Ui are the baryon density, the pressure, the specific energy, the three velocity, and the spatial components of the 
four velocity, e is the sum of the specific internal energy e and the specific kinetic energy as 

e = e+^UiUi. (4) 

and Ui are related by 

Ui=lijV^, (5) 
where 'jij = Sij + hij . For consistency, we compute by 

yi=f3Uj, (6) 

where 7*'' is the inverse of 7^. 

tl), ipR, and hij denote the Newtonian potential, the radiation reaction potential, and the tracefree tensor component 
of the radiation reaction metric. V and V'fl are determined from the equations 



Alp = 477/9, (7) 



fc^lf-J + Aar-Ul. (8) 



where $ obeys 

A di = Am-h ■ ■ T-J - _ 

dx^ 



A^=47rhijx^^. (9) 



hij can be related to the tracefree quadrupole moment -l-ij as 

4G d^4- 



where 

^ij = j d^xp(^x'x^ - ^r^S'^y (11) 

In Eq. (10), we recover c to clarify that it appears only in hij-. In other terms of the basic equations, c is not 

included explicitly, e is introduced to control the strength of the radiation reaction force. Note that tjjji and <I> are also 
proportional to e. Hence, the radiation reaction time scale will be proportional to e~^. To confirm that the growth 
time scale of the secular instability driven by gravitational radiation is indeed proportional to the inverse of e, we 
performed simulations varying the value of e. 

The order of magnitude for d'' -I- ij / dt^ for rapidly rotating stars is estimated as ~ (GM j Rf?)^!"^ . Thus, hij is totally 
proportional to e{GM / Rc'^)^/'^ . This implies that increasing the value of e is equivalent to increasing the value of the 
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compactness in the (0+2.5) post-Newtonian formalism. We note that the compactness in this framework does not 
have any meaning with regard to the strength of the self-gravity since in this framework, c appears only in the terms 
associated with the gravitational radiation reaction. The compactness here is only the indicator for the magnitude of 
the radiation reaction force. 

The third time derivatives of -J-ij cannot be computed accurately by the straightforward finite- differencing since we 
adopt a second-order accurate finite-differencing in time. For the accurate computation, we adopt the same method 
as that adopted in previous papers (e.g., [37,38]). In this method, the time derivatives of the quadrupole moment are 
rewritten appropriately operating the time derivatives in the integral and using the equations of motion. Then, the 



third time derivative of rf ^ is written 



i_^ = STF 



ox' ox^ ox3 ax" ox^ 



0(e), (12) 



where STF implies that the symmetric tracefree part should be extracted as 



STF[T,,] = ^ii±^-^Tfe,, (13) 

and Tij is a tensor. The terms of order e in d^-I-ij/dt^ yield the terms of order in hij, and thus, it is neglected. 
During numerical simulation, we adopt the F-law equation of state 

P = (F - l)pe, (14) 

where F is the adiabatic index, and in this paper, we set F = 2 to model a moderately stiff equation of state for 
neutron stars as in [17]. With this equation of state, we can follow the formation of shocks that may form after the 
secular instability grows to a nonlinear regime. 

The hydrodynamic equations are solved using a so-called high-resolution shock-capturing scheme we have recently 
adopted in fully general relativistic simulations [39,35]. In this method, the transport terms in the hydrodynamic 
equations di{- ■ •) arc computed by applying an approximate Ricmann solver with third-order (picccwise parabolic) 
spatial interpolation. Detailed numerical methods with respect to the treatment of the transport terms (in the 
framework of general relativity) are described in [39]. The time integration is done with the second-order Runge- 
Kutta method. Atmosphere of small density 10~^ of the central density pc) outside stars is added for using the 
shock-capturing scheme. Since the hydrodynamic equations are written in the conservative form, the mass is conserved 
within the digit of double precision. The conservation of energy and angular momentum is slightly violated due to 
numerical dissipation besides the loss by gravitational waves. 

In the (0-1-2.5) post-Newtonian formulation adopted here, we need to solve three Poisson equations for ^, dti^, 
and They are solved by a preconditioned conjugate gradient method, for which the details are given in [38], with 
boundary conditions 

V' - - ^^ijnV + 0{r-% (15) 

5tV ^ - ^3 + Oir-'), (16) 

^^-^ajn^n^ + 0{r-% (17) 

where n' is a radial unit vector x'/r, Cij is a moment defined from the right-hand side of Eq. (9), and M is the total 
baryon mass defined by 



M 



= J d^xp. (18) 



We adopt a fixed uniform grid with size 141 x 141 x 141 in x — y — z, which covers a region —L < x < L, —L < y < L, 
and < z < L where L is a constant. We assume a reflection symmetry with respect to the equatorial plane, and set 
the grid spacing of z to be half of that of x and y. In all the simulations, the equatorial radius of rotating stars 
is set to be 5L/7 (i.e., the equatorial radius is covered by 50 grid points). We also performed test simulations with 
size 113 X 113 x 113 (i.e., the grid spacing becomes 1.25 times larger) for several selected cases and confirmed that 
the results depend weakly on the grid resolution. In the absence of radiation reaction, we checked the conservation 
of the total energy and angular momentum. The magnitude of the error monotonically increases with the number of 
time steps and for a typical simulation with 141 x 141 x 141 grid size, after 30 central rotational periods {Pq), the 
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violation of the total energy and angular momentum conservation is ^ 5% and 2%, respectively. We have checked that 
these values are decreased at approximately second order with improving the grid resolution. Even in the presence 
of gravitational radiation reaction, the vorticity should be conserved. We have monitored its conservation and found 
that for a region of relatively high density ^ 10~^Pc where pc denotes the central density, it is conserved within 
~ 10% error. (Around the surface of stars, the error is larger.) The violation of these conservation relations is due to 
numerical dissipation or diffusion. 

Numerical simulations were performed on FACOM VPP5000 in the data processing center of National Astronomical 
Observatory of Japan. For one numerical simulation with size 141 x 141 x 141, it takes about 8 CPU hours for 80,000 
time steps using 48 processors. 



B. Initial conditions 



A slightly perturbed rotating star is adopted as the initial condition for numerical simulation. Specifically, we first 
construct models of differentially rotating stars in axisymmctric equilibrium with values of /? = T/|M^| past the secular 
instability threshold [12]. Here, in the notation of this paper, T and W are computed from 



T 

W = 



^ j (fxfruj'^n'^, (19) 
= \j Sxfnp. (20) 



The equilibrium state is computed using the polytropic equation of state, 

P = Kp^, (21) 

where K is the polytropic constant and F = 1 + 1/n with n polytropic index. In this paper, we choose F = 2 (n = 1). 

Since a realistic velocity profile for protoneutron stars is poorly known (e.g., [23]), we give a plausible one. As the 
angular velocity profile Q{w) where w = \/ x"^ + j/^, we choose the so-called j-constant-like law as 

where A is a constant, and the angular velocity at the rotational axis {z axis). The parameter A controls the 
steepness of the angular velocity profile: For smaller values of A, the profile is steeper and for A ^ oo, the rigid 
rotation is recovered. In the present work, the value of A is chosen to be equal to the equatorial radius R,,. In this 
case, the ratio of f2o to Q at the equatorial surface is two, and hence, the star is moderately differentially rotating. 
Numerical simulations have shown that such a rapidly rotating star of moderate degree of differential rotation can 
be formed from initial conditions with rapid and moderately differential rotation [23] . As we have shown in previous 
papers [30], a star of a high degree of differential rotation, i.e., A = A/Re ^ 0.5 is dynamically unstable against 
nonaxisymmetric deformation even if it is not rapidly rotating with /? = 0(0.01). However, with A = 1, the star is 
dynamically stable for /3 ^ 0.265 [30,40], and hence, the rapidly rotating stars with 0.2 < j3 < 0.26 that we adopt in 
this paper (cf. Table I) can be subject only to the secular instability. We also note that only for A ^ 1, a rotating star 
of a high value of /? ~ 0.25 can be found for F = 2 (e.g., [12]). Such a star which is secularly unstable but dynamically 
stable may be a plausible initial state of a rapidly rotating protoneutron star as mentioned in Sec. I. In addition, only 
for such a high value of /3, the growth time scale of a secularly unstable mode is short enough to accurately follow the 
nonlinear growth of the perturbation. 

In terms of /? and A, one rotating star is determined for a given rotational profile and polytropic index. To specify 
a particular model, we could also specify the ratio of the polar radius Rp to the equatorial radius Re, i.e. Ca = Rp/Rg 
instead of /3. For the equations of state and the angular velocity profiles that we adopt in this paper, the value of Ca 
monotonically decreases with increasing /J for any given set of A and n. This is the reason that Ca can be a substitute 
for /3. 

We initially add a secularly unstable perturbation of a small amplitude on top of the equilibrium states. The 
density and velocity profiles of an unstable mode adopted are those computed by Karino and Eriguchi [12] in a linear 
perturbation analysis. In their analysis, the perturbed quantities are expanded as 

5p = Spmjr, 0) cos{m,(p — at), (23) 
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(a) (b) 

FIG. 1. (a) The profiles of the density, the rotational velocity, and the perturbed quantities on the equatorial plane, and 
(b) the density contour curves in x-z plane for Ca = 0.3. The density and the rotational velocity are normalized by pc and 
Q()Rc, respectively. The amplitude of the perturbation is normalized so that the maximum value of is unity in units of 
G = pc = Re = Here, v"^ = wv^ . The density contour curves are drawn for p/pc = O.lj for j = 0.01, 0.1, 1, 2, • • • , 8, 9, 9.5, 
and 10. 

Sv"" = ^ Sv^ (r, 9) sm{mip - at) , (24) 

m 

6v^ = ^5vi,{r,9)siii{mip-at), (25) 

m 

Sv"^ = ^ 6vg{r, 6) cos{mip - at), (26) 



where a is the angular velocity of the perturbation. The eigenfunctions are these found by solving the perturbed 
equations for Spm and Sv}^ in the two dimensional space of r and 6. Here, and v'^ denote the components associated 
with the orthonormal bases. We put an unstable mode of m = 2 on a rotating equilibrium star at t = 0. The profile 
of the perturbed quantities as well as the density and the rotational velocity on the equatorial plane are shown in Fig. 
1 for Ca = 0.3 for m = 2, the only mode considered here. Note that the profile of the perturbed quantities is not as 
monotonic as that for rigidly rotating and incompressible stars [1]. Figure 1(b) displays density contour curves in x-z 
plane for Ca = 0.3 and shows that the rotating star is a flattened spheroid. In Table I, we also list several quantities 
of the perturbed rotating stars adopted as the initial condition. Here, the energy and the angular momentum are 
defined by 

E = j d^xpe + W, (27) 
J = j (fxpu^. (28) 

The growth time of a secularly unstable mode is approximately proportional to e^^ = er^ {GM / ReC^)^^^"^ ■ Thus, 
to accelerate the growth, we artificially increased the value of e by a factor of 10-100. The values chosen in this work 
are listed in Table L In the following, we will refer to the case listed in Table I as "e = 1" case. To increase e, we may 
either increase the value of e or the value of compactness. As mentioned in Sec. II A, the compactness can be rather 
arbitrarily increased in the (0+2.5) post-Newtonian framework. In our choice of e for numerical computation, setting 
e = 1 is equivalent to choosing very compact stars with GM/RgC^ 0.6, which is about 3-6 times as large as that of a 
neutron star or protoneutron star. We note that in general relativity, a star with GM/RgC? ^ 0.6 collapses to a black 
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T/\W\ 


M/{pMt) 


E/iM^R^) 


j/{m^/^rY') 




cr/f2() 




I^^/{MRl) 


e 


r/Po 


Tanal/ Po 


0.30 


0.256 


0.567 


-0.538 


0.366 


1.34 


0.569 


0.914 


0.288 


0.25 


3.5 


5.9 


0.334 


0.237 


0.557 


-0.546 


0.347 


1.30 


0.455 


1.02 


0.278 


0.24 


15 


34 


0.35 


0.229 


0.565 


-0.547 


0.339 


1.29 


0.419 


1.06 


0.275 


0.25 


31 


60 


0.40 


0.203 


0.610 


-0.547 


0.315 


1.26 


0.290 


1.16 


0.272 


0.31 




510 



TABLE I. Several quantities for secularly unstable rotating stars are shown in units oi G = Re = pc = 1 where pc is the 
central density at t = 0. cr is the angular velocity of the unstable mode and r is the growth time scale of the secularly unstable 
mode found in numerical simulation for the value of e = e{GM / Re.(?Y'^'^ given here in units of Po = 2-k/Q.(, (the rotational 
period along the z aocis). Note that the value of e of a realistic neutron star with the equatorial radius 10-20 km and the mass 
1.4M© is 0.018-0.003. " — " implies that we could not determine the value of r since it is too large. For all the models, the 
central density is normalized to be unity. 



hole. However, in the (0+2.5) post-Newtonian theory, a star of any value of GM/R^c? is stable against collapse for 

r = 2. As a result, we can adopt such a large compactness of which the growth time scale of the secularly imstablc 
mode is short enough to accurately follow its growth numerically. This is a great advantage in this post-Newtonian 
framework. 



C. Method for analysis 

The growth of a nonaxisymmetric bar-mode perturbation is monitored using a distortion parameter defined as 

^^{rjl + r^lf'\ (29) 

where 

ri+ = (30) 

J-xx ~r J-yy 

and lij = x,y, z) denotes the quadrupole moment. 

Throughout this paper, we choose initial conditions in which r] ~ 0.01 or 0.05 at f = by appropriately multiplying 
a constant to 5pm and (^u^. (Hereafter, we denote the initial value of rj as 779.) When the instability grows, 77 should 
increase in proportional to e*/"^ where r denotes the growth time scale. Thus, we fit the time evolution of rj and 
determine the value of r. 

During numerical simulation, gravitational waves are computed in the quadrupole formula [41] in which "+" and 
"x" polarization modes of the waveforms are defined as 

h+ = bl^Im^ = (32) 

r r 

and the luminosity and the angular momentum flux by 

^w^^E(^S^S^-^S^^S)> (34) 

3 

where 

4 = ^, ^.(f = ^ (n = 2,3), (35) 



and r is the distance from the source to a detector. h+ and hx presented here are the waveforms observed along the 
rotational [z-) axis. Thus, eff'ectively, they represent the maximum amplitude. 
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(a) (b) 
FIG. 2. (a) as a function of time for Ca = 0.3 (solid curve), 0.334 (dashed curve), and 0.35 (dotted-dashed curve) with 
€ = 1 and rja ~ 0.01. For comparison, the results with rja ~ 0.04 are shown for Ca ~ 0.334 (long-dashed curve) and Ca = 0.35 
(dotted-long-dashed curve). The dotted curves denote the results with a low-resolution grid for Ca ~ 0.30 and Ca ~ 0.334 with 
rjo ~ 0.04. (b) r; as a function of time is shown for Ca = 0.3 with e = 1 (solid curve), 0.7 (long-dashed curve), 0.5 (dotted-dashed 
curve), and (no radiation reaction; dotted-long-dashed curve). The dotted line denotes rj = rjo. For e = 0.5 and 0.7, results 
with rjo ~ 0.05 are shown together. 

III. NUMERICAL RESULTS 

A. The growth of secularly unstable bar-mode perturbations 

Numerical simulations were performed for secularly unstable stars with Ca = 0.3, 0.334, 0.35, and 0.4. For the 
case with Ca ^ 0.52 (/3 ^ 0.14), the star of A = 1 and F = 2 is secularly unstable according to the result of a linear 
perturbation analysis by Karino and Eriguchi [12]. In Fig. 2 (a), we show 77 as a function of time for Ca — 0.3, 0.334, 
and 0.35 with e = 1. For all cases, rj increases with time approximately as e*/'^. The plausible reason for a slight 
deviation from the exponential form is that modes other than the secularly unstable ones are excited due to numerical 
noise or nonlinear coupling, which contribute to the growth and modulation of 77 soon after we start the simulation. 
However, they are not unstable modes, and hence, their amplitude does not increase. As a result, after the amplitude 
of the unstable mode increases for t ^ 2Po, the effect of other modes does not seem to be very large. From this reason, 
we determine the value of r from the curve of 77 for i ^ 2Po- More specifically, we use the data sets of 2Po <t< TPq 
for Ca = 0.3, of 2Po<t< 20Po for Ca = 0.334 with rjo ~ 0.04, and of 2Po < t < 30Po for Ca = 0.35 with rjo - 0.04. 
We expect that the error for the evaluation of r is ~ 10% for Ca = 0.3 and 0.334. For Ca = 0.35, the growth time scale 
is too long to follow the late evolution with a sufficient accuracy. Moreover, a modulation in ij prevents the accurate 
evaluation of the growth time scale. Thus, the magnitude of the error in the numerical result of t for Ca — 0.35 may 
be > 10%. 

From the analysis, we determine t/Pq ~ 3.5, 15, and 31 for Ca ~ 0.30, 0.334, and 0.35 (see Table I). For Ca = 0.334 
and 0.35, the growth time is also evaluated from the data with 770 ~ 0.01, and we find that the results agree 
approximately with those for 770 ~ 0.04. This shows that the growth time is independent of the value of t^q. 

According to a linear perturbative analysis for the bar-mode secular instability of incompressible and rigidly rotating 
stars, the growth time of the secularly unstable perturbation induced by emission of gravitational waves is [5] 

T-anal = , (36) 

where a is the angular velocity of the secularly unstable mode, (T+ is that of its conjugate mode, which is a secularly 
stable one, and /^ra = Ixx + iyy Equation (36) is valid only for the incompressible and rigidly rotating case. For 
7T, ^ 0, it may be an approximate formula [4], but it is not very clear. Using the numerical results for a and cr+ in 
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a linear pcrturbativc analysis by Karino and Eriguchi [12], the value of Tanai is evaluated and listed in Table I. It is 
found that the value of t computed from the numerical results is systematically smaller than Tanai^ '''/''anal ^ 0.5 0.6. 
The systematic disagreement will be due to the fact that Eq. (36) is valid only for the incompressible and rigidly 
rotating star. Indeed, the profile of the perturbed quantities is significantly different from that for incompressible and 
rigidly rotating stars [1]. On the other hand, the value of r increases steeply as a/flo decreases by a similar manner to 
that described in Eq. (36). This indicates that the perturbation certainly grows due to the secular instability induced 
by gravitational radiation reaction. 

To check that r is proportional to e~^, in Fig. 2 (b), we show 77 as a function of time for Co = 0.3 with e = 1 (solid 
curve), 0.7 (long-dashed curve). 0.5 (dotted-dashed curve), and (no radiation reaction; dotted-long-dashed curve). 
In these simulations, the compactness GM/RgC^ is fixed. For e = 1, % ~ 0.01, but for e = 0.7 and 0.5, ryo is set to be 
~ 0.01 and 0.05. For the case with no radiation reaction, the value of 77 does not increase but simply oscillates due 
to the growth of modes other than the secularly unstable one probably by numerical noises and nonlinear coupling. 
This shows that the star with Ca — 0.3 is dynamically stable. (For A = 1 and F = 2, only the stars with Ca ^ 0.25 
are dynamically unstable as demonstrated in [30].) The growth time r is determined by the same method as that 
adopted above. The analysis gives t/Pq = 3.5, 4.9, and 7.4 with an error of magnitude ~ 10% for e = 1, 0.7, and 
0.5. Thus, t/Pq is approximately proportional to e~^. We also note that r also depends very weakly on the value of 
rjo <; 0.05. This implies that the amplitude of the nonaxisymmetric perturbation initially given is small enough. 

Figure 2(b) also indicates that the value of r] eventually reaches ~ 0.3 and remains between ~ 0.2 and ~ 0.3. This 
shows that an ellipsoid is formed. The maximum value of 77 (hereafter rymax) at the formation of the ellipsoid depends 
weakly on the magnitude of the parameter e. However, the shape of the curve of rj{t) for rj ^ 0.1 is fairly different. 
For e = 1, r] increases almost monotonically to the maximum value, but for e = 0.7 and 0.5 with 770 ~ 0.01, the growth 
rate is varied at a time that 77 becomes ~ 0.1. A plausible explanation for this difference is that the radiation reaction 
time scale is too short to obtain a qualitatively independent result of the value of e: We recall that the radiation 
reaction force adopted here has an invariant meaning only when the system is adiabatic, i.e., the system is nearly 
periodic and the reaction time scale is much longer than the period of a nonaxisymmetric oscillation. For Ca = 0.3 
with the compactness ~ 0.6, the reaction time scale is only a few times as long as the period of the nonaxisymmetric 
oscillation ~ 2Po ~ Te~^/2. This fact may be reflected in the difference of the shape of the curve of 77. 

In Fig. 2 (a), we display the results with a low-resolution grid for Ca = 0.3 and 770 ~ 0.01 and for Ca = 0.334 
and ?7o ^ 0.04 (dotted curves). For both cases, e = 1. The growth time of 77 and the maximum value of 77 agree 
approximately with those for the high-resolution simulation. The growth time is slightly longer for the lower resolution 
simulations. This indicates that for a larger numerical dissipation, the growth time scale is spuriously increased. 
However, the grid resolution adopted here is high enough to obtain a convergent result. 

B. The fate of unstable stars 

In Figs. 3 and 4, we display snapshots of the density contour curves and the velocity vectors in the equatorial plane 
for Ca = 0.3, e = 1, and 7/0 « 0.01, and for Ca = 0.334, e = 1, and 770 ~ 0.04. With the growth of a secularly unstable 
and nonaxisymmetric perturbation, the spheroidal shape is changed to an ellipsoidal shape, and after the value of 77 
reaches 0.1 0.3, the nonlinear growth of the perturbation terminates. Since the value of r\ remains of order 0.1, the 
termination of the growth of the nonaxisymmetric deformation appears to happen when the rotating star reaches a 
quasistationary ellipsoidal state (see below for more discussion). Here, "quasistationary" implies that the axial ratio 
of the ellipsoid does not change in a dynamical time scale and the evolution is determined by the radiation reaction 
time scale which is much longer than the dynamical time scale of system. 

The value of r/max does not depend strongly on the value of 770 [cf. Fig. 2(a)], but it docs on Ca- In particular, the 
smaller the value of Cq, the larger the value of 77max so that 77max ~ 0.3, 0.2, and 0.1 for Ca = 0.3, 0.334, and 0.35, 
respectively. This feature is expected from a study for the incompressible fluid [16] and as well as from a compressible 
ellipsoidal model [4] . As the maximum value of r\ is reached, a semi major axial ratio of the formed ellipsoid defined in 
the equatorial plane also depends on the value of Ca- For Ca = 0.3, the axial ratio of the ellipsoid ~ ^(1 — 77)/(l + 77) 
reaches ~ 0.75, and the numerical results suggest that it is smaller for the larger value of Ca- We note that for 
an incompressible star [4], the axial ratio reached for /? ^ 0.25 should be much smaller than 0.75. This may be a 
consequence that a rotating star with larger value of n is less prone to forming an ellipsoid [42]. However, even for 
F = 2, the axial ratio of the semi major axes is fairly large for Ca ^ 0.35, and therefore, the ellipsoid will subsequently 
emit gravitational waves of large amplitude and dissipate energy and angular momentum due to a longterm emission 
of gravitational waves (cf. Sec. Ill C). 

In Fig. 5, we show the profiles of the density and angular velocity along the x and y axes at selected time slices. 
The contour curves and the velocity vectors at the corresponding time steps are shown in the third panels of Figs. 3 
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and 4. In other words, Figures 5(a) refers to the profile at t/Po = 12.8 for Ca = 0.3 and e = 1, and Fig. 5(b) is at 
t/Po = 24.9 for Ca = 0.334 and e = 1 with rjo w 0.04. In each figure, the density and the angular velocity at t/Po — 
(dotted curves) are shown together. Around the origin, the density profile does not change significantly even after the 
formation of the ellipsoid. (Note that the central density decreases with time partly due to a numerical dissipation 
or diffusion.) However, around va/Re ~ 0.5-0.6, a density peak along the semi major axis is formed irrespective of 
the initial condition. Also, it is found that a nonaxisymmetry is enhanced around zu/Re > 0.2. These features found 
in the density profiles reflect the initial perturbation profile [cf. Fig. 1(a)]. After the formation of the ellipsoid, 
differential rotation is also enhanced as the result of the growth of the nonaxisymmetric density perturbation. It is 
interesting to note that a high-density and nonaxisymmetric peak rotates with a smaller rotational speed than that 
for the unperturbed axisymmetric star. This seems to be reasonable because the secularly unstable mode has an 
angular velocity smaller than the rotational angular velocity of the unperturbed equilibrium state. 

In Fig. 6, the evolution of the total angular momentum J and P are shown for Ca = 0.3, and e = 1 (solid curve), 
0.7 (long dashed curve), and (dashed curve). Even for e = for which the system should be almost stationary 
throughout the simulation, J increases by ^ 2% and f3 decreases by ~ 1% in 30Po. These are the typical truncation 
error with grid of 141 x 141 x 141. For e > 0, during the exponential growth of the secularly unstable mode, J and 
/? decrease at most 1% by gravitational wave emission. (Compare the curves with e = and e ^ 0. The difference 
between two values can be regarded as loss by the gravitational radiation.) Thus, the formed ellipsoid initially has 
slightly smaller values of J and (3 than those of the secularly unstable axisymmetric star in equilibrium. After the 
formation of a quasistationary ellipsoid, these values start decreasing significantly due to dissipation by gravitational 
waves. Figure 6 indicates that when the ellipsoid relaxes to a final stationary state that does not emit gravitational 
waves, the values of J and /3 will be much smaller than the initial values. 

The final state after a sufficient emission of gravitational waves is uncertain. For an incompressible fluid star, it will 
be a Dedekind ellipsoid which is static in shape and has only internal motion and no pattern rotation if it is observed 
in the inertial frame. However, for n 0.8, such configuration may not exist because of the following facts: (i) the 
uniformly rotating ellipsoid (Jacobi-like ellipsoid) does not exist for n > 0.808 [42] and (ii) the Dedekind theorem for 
the incompressible fluid tells that the absence of the Jacobi ellipsoid implies the absence of the Dedekind ellipsoid [1] . 
Namely, if this theorem holds for n 7^ 0, no Dedekind-like configuration would exist for n > 0.808. Lai and Shapiro 
expect that the final state may be a Dedekind-like ellipsoid even for n > 0.8 [4]. The other possibility is that the 
unstable star does not relax to a Dedekind-like ellipsoid but to a (nearly) axisymmetric and secularly stable rotating 
star after sufficient emission of gravitational waves. The results of the present simulation indicate that the latter is 
more likely because of the following reasons. 

For the incompressible case, a secularly unstable Maclaurin spheroid changes to a Dedekind ellipsoid by emission of 
gravitational waves conserving the vorticity [4]. Throughout this transition, the energy and the angular momentum 
are decreased at most by a few perccnts from their initial values [4]. Furthermore, for P ~ 0.25, the axial ratio in 
the equatorial plane becomes very small ^ 0.75. In the case of Ca — 0.3, 770 ~ 0.01, and e = 1, until the end of the 
simulation at t ^ 30Po) the energy and the angular momentum are dissipated by gravitational waves by ^ 5% and 
~ 10%, respectively. However, even at f ~ 30Po, the rotating stars still do not settle into a highly deformed Dedekind- 
like ellipsoid but to a moderately deformed quasistationary (Riemann-type) ellipsoid which emits gravitational waves 
of a moderately large amplitude (cf. Sec. Ill C). The results of the simulation also indic;ate that the gravitational 
wave emission does not enforce the ellipsoid to a stationary state immediately, although the final outcome should 
be an object that does not emit gravitational waves. Thus, the ellipsoid will subsequently emit gravitational waves 
for a long time to dissipate the energy and the angular momentum by <; 10%. Indeed, Fig. 6 suggests that f3 will 
decrease much below the initial value. Therefore, in the present case, the final state after a longterm dissipation 
by gravitational waves will not be a highly deformed Dedekind-like ellipsoid, but a nearly axisymmetric spheroid of 
/3< 0.25. 

As mentioned above, the uncertainty on the existence of the Dedekind ellipsoid for compressible stars may be a 
key issue in understanding the final state after a longterm emission of gravitational waves. To investigate whether 
a compressible Dedekind equilibrium exists or not, a specially designed implementation is necessary. Unfortunately, 
systematic numerical computation of the compressible Dedekind ellipsoid has not been succeeded yet because of the 
technical difficulty (but see [43]). A robust numerical implementation is required to clarify the possible final state 
after a longterm dissipation by gravitational waves. 

C. Gravitational waves 

In Fig. 7, the gravitational waveforms, luminosity, and angular momentum flux for Ca = 0.3 with e = 1 are shown. 
In the early stage in which rj increases exponentially, the amplitude and the fluxes of gravitational waves are increased 
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in an exponential manner. After the end of the nonlinear growth of the nonaxisymmetric perturbation, they relax 
approximately to constant values. (The reason that they do not exactly settle into constants is that modes other than 
the secularly unstable mode are excited.) This global feature of gravitational waveforms is independent of the initial 
value of Ca, although the amplitude is different depending on the value of 77. 

When the secularly unstable mode of the angular velocity a is the only one present mode, the luminosity of 
gravitational waves should be equal to 

iGwi = ^<yWil^- (37) 

To confirm that the amplitude of gravitational waves is increased mainly by the growth of the secularly unstable mode, 
we compare the energy flux -Lqwi with i^Gw- In Fig- 7(b), the evolution of Lgwi is shown together (dashed curve). 
We see that the luminosities computed by two methods agree approximately, although in the curve corresponding to 
Lgw shows oscillations. This is due to the fact that modes other than the secularly unstable one are excited. 

From a break of the luminosity curve at t ~ 12Po) we can approximately identify the end-point of the exponentially 
growing stage of the secularly unstable perturbation. Until this time, the energy and the angular momentum carried 
away by gravitational waves are only by ^ 1% of the total energy and angular momentum for all the models. This 
implies that the quasistationary ellipsoid formed at the end of the growth of the perturbation has only slightly smaller 
energy and angular momentum than those of the secularly unstable and axisymmetric rotating star. In contrast to 
the incompressible case [4], most of the rotational kinetic energy and the angular momentum are dissipated by a 
subsequent emission of gravitational waves from a formed ellipsoid as discussed in Sec. Ill B. In the following, let us 
estimate the accumulated gravitational wave amplitude assuming that the ellipsoid is a protoneutron star. 

For the luminosity Lew = 0.001a(M/i?e)^ and the kinetic energy T = ({M^/Rg) where a and ( are parameters, 
the emission time scale of gravitational waves of a protoneutron star can be estimated as 

12«-fAV:^VfT^V'3ec. (38) 



_0.2yV20kmy \1AMq_ 

Here, a w (t//0.3)^ ^ 1 and ( ^ 0.2 for (3 = 0.2 0.25, respectively. The value of Re will depend strongly on 
the rotational profile of a progenitor of supernova stellar collapse (e.g., [20]). In the present case, we assume that 
protoneutron stars are rapidly rotating with /3 > 0.2. In such case, the equatorial radius would be larger than a 
canonical radius of neutron stars ^ 10 15 km because of the strong centrifugal force. For a sufficiently large initial 
value of /? ^ 0.01 and for a small value of A <, 1/2, the outcome of stellar collapse could be an oscillating star of 
subnuclear density and of a high value of /3. In such case, Re may be as large as ~ 100 km. 
The characteristic frequency of gravitational waves is denoted as 

1 /2 

where a is the angular velocity of the nonaxisymmetric perturbation in units of pc , where pc is the initial central 

1 /2 

density, cr/pc depends strongly on the value of /3, but the value of a is always smaller than by a factor of ^ 2. 
This is a very important feature for the detection of gravitational waves by kilometer-size interferometers for which a 
sensitivity is the best in the frequency band around a few 100 Hz <C i^o/i"- 

Assuming that the nonaxisymmetric perturbation would not be dissipated by viscosities or magnetic fields on the 
emission time scale of gravitational waves [44] , the accumulated cycles of gravitational wave-train N are estimated as 

In reality, the characteristic frequency will be changed as gravitational waves are emitted, and hence, the expression 
of N which denotes the accumulated cycle for a given frequency should be regarded as the maximum value. However, 
the order of magnitude for the real value of A'' will be identical. 

The effective amplitude of gravitational waves is defined by hcg = N^^'^h [45,4,46] where h denotes the characteristic 
amplitude of periodic gravitational waves. According to the numerical results. 
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where h ocr] and ~ a^^^ ^ 1. Using this relation, 

Here, a~^/'^h ^ 1 irrespective of 77, and thus, /icfr is approximately independent of 77. This property is qualitatively 
explained as follows: As gravitational waves are emitted, the value of 77 will be decreased. This implies that the 
amplitude of gravitational waves h and luminosity are decreased. However, the time scale of gravitational radiation 
reaction is increased with decreasing the luminosity. These two effects cancel each other as that /icfr is approximately 
independent of 77. Thus, Eq. (42) may be used for any value of 77, i.e., for any initial value of /?. (Remember that the 
maximum value of 77 is smaller for smaller values of 

Equation (42) shows that h^s will be larger than IQ-^^ at a distance r 50 Mpc with the frequency of about 500 
Hz for Re ~ 20 km and M k. IAMq. As gravitational waves are emitted, the rotational velocity of the ellipsoid will 
be decreased, and hence, the frequency of gravitational waves and the value of C, will be also decreased. However, 
hes is proportional to ct^/^ and C^^^ for given values of mass and radius, and thus, the dependence of h^s on these 
two parameters is not very strong [4]. Therefore, the value of /leg could be ^ IQ-^^ at the frequency at a distance 
of ^ 10 Mpc in a wide range between ^ 10 Hz and ^ 500 Hz for ^ 20 km and M « 1.4M0. (This amplitude 
agrees approximately with that derived in [4].) This implies that gravitational waves from protoneutron stars of a 
high value of /3 and of mass ~ 1.4M0 and radius ~ 20 km at a distance of ~ 10 Mpc are likely to be broadband and 
quasipcriodic sources for laser interferometric detectors such as LIGO [47]. 

If other dissipation or transport processes of the angular momentum, e.g., by viscosity or magnetic field, is present, 
growth of the nonaxisymmetric perturbation in a protoneutron star may be suppressed. A detailed discussion about 
the viscous effects is found in [4], in which the authors indicate that the viscous effect will be negligible. On the 
other hand, effects due to magnetic fields have not been investigated in detail. In [44], the authors estimate the 
time scale of magnetic braking which transports the angular momentum in a differentially rotating star outward in 
~ 10(i?/10^2Gauss) sec where i? is a strength of the magnetic field. For a canonical value oi B = 10^^ Gauss, the 
time scale is comparable to Tevo- Therefore, in the presence of the magnetic fields, a rotating ellipsoidal star may not 
evolve simply due to gravitational radiation reaction. This point is not clear at present. 

IV. SUMMARY 

We have numerically studied the growth of a secularly unstable bar-mode induced by gravitational radiation for 
rapidly rotating protoneutron stars of /? between ~ 0.2 and ~ 0.25 with the F = 2 polytropic equation of state. 
Numerical simulations were performed in a (0+2.5) post-Newtonian framework in which the effect of gravitational 
radiation reaction is incorporated in addition to Newtonian gravity. In our simulations, the amplitude of the secularly 
unstable bar-mode is driven to a nonlinear stage in which the deformation parameter becomes 0.1-0.3 depending 
on the initial value of /3. The outcome is a moderately deformed ellipsoid of the semi major axial ratio ^ 0.75, in 
contrast to the incompressible case, in which it becomes much smaller than 0.75 for [3 ~ 0.25. The final values of 
the deformation parameter 77 and of the axial ratio do not depend strongly on the radiation reaction force. However, 
they depend on the initial value of /?, and are smaller for the smaller value of /?, as in the incompressible and rigidly 
rotating case [16,4]. The growth of the nonlinear perturbation appears to be terminated because the formed ellipsoid 
is quasistationary. The total energy and angular momentum dissipated by gravitational waves until the formation of 
the ellipsoid is ~ 1% of the initial values. Thus, the formed ellipsoid initially has only slightly smaller energy and 
angular momentum than those for the axisymmetric rotating star adopted at t = 0. 

In contrast to the incompressible and rigidly rotating case, the formed quasistationary ellipsoid is not likely to 
be a Dedekind-like stationary ellipsoid. Namely, their pattern rotates globally, and hence, gravitational waves are 
emitted. The simulations indicate that the gravitational wave emission does not enforce the ellipsoid to a stationary 
state immediately, although the final state should be such a state that does not emit gravitational waves. Thus, 
the ellipsoid will subsequently emit gravitational waves for a long time to dissipate the energy and the angular 
momentum significantly. This implies that the final state may not be a highly deformed Dedekind-like ellipsoid but 
a nearly axisymmetric spheroid. Also, this indicates that in the absence of other dissipative mechanisms, a rapidly 
rotating protoneutron star may be a stronger emitter than that considered before [4], because the energy and the 
angular momentum will be dissipated by gravitational waves by ^ 10%, not a few %, until formation of a stationary 
star. 
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The growth time of rj (denoted by r) is cvahiated using numerical data sets and compared with an analytic formula 
Tanai for incompressible and rigidly rotating stars. It is found that the value of t is systematically smaller than 
that derived from Tanai- Plausible explanations for this difference are either of the following facts: (i) the analytic 
formula is valid only for rigidly rotating stars and, hence, in differentially rotating stars, r becomes systematically 
shorter than Tanai; (h) the analytical formula is valid only for incompressible stars and for compressible stars, r 
becomes systematically shorter than Tanai- To clarify the plausible reason, a linear perturbative analysis including 
gravitational radiation reaction will be the best approach. 

An ellipsoidal neutron star formed after the growth of a secularly unstable mode of a rapidly rotating progenitor 
has a fairly large ellipticity and a pattern rotation. The angular velocity of the pattern rotation is much smaller 
than that of a axisymmetrically rotating equilibrium star initially given. Thus, the formed ellipsoid will be a strong 
emitter of gravitational waves of fairly low frequency, less than several hundreds Hz. Using a simple analysis, an 
effective amplitude of gravitational waves is evaluated. It is found that for a protoneutron star of the equatorial 
radius Re ~ 20 km, mass M w 1.4Mq and the initial value of /3 ~ 0.25, the effective amplitude hctf will be larger 
than 10~^^ at a distance r ~ 10 Mpc in a wide frequency range between 10 and ^ 500 Hz. Such gravitational 
waves from protoneutron stars are likely to be sources for laser interferometric detectors such as LIGO [47]. In the 
presence of magnetic fields of magnitude ^ 10^^ Gauss, magnetic braking may affect the evolution of differentially 
rotating and ellipsoidal protoneutron stars significantly. In such case, the evolution for the ellipsoidal protoneutron 
star is not yet clear at present, and thus, a further study is necessary. 

Note added in proof. Soon after this paper was submitted, a paper by Ou, Tohline, and Lindblom [48], which also 
studies the secular instability driven by gravitational radiation, was posted in astro-ph. They performed a (0+2.5) 
post-Newtonian simulation for a rapidly rotating polytropic star with F = 3 including a radiation reaction term in 
essentially the same approximate manner as that adopted in [17]. They find the similar result to that found in our 
present paper for the growth of the secular instability to form an ellipsoid. The formed ellipsoid has a pattern rotation 
of small velocity, and thus, is not the Dedekind but the slowly rotating Riemann-type ellipsoid. On the other hand, 
they chose F = 3 that is larger than our choice (F = 2). As a reasonable result, they find that the axial ratio of the 
formed ellipsoid is slightly smaller (~ 0.5) than that found in our work. The radiation reaction formalism we adopt 
can be used for any problem as long as the system adiabatically evolves due to the back reaction. In contrast, the 
radiation reaction formalism they adopt can be in principle used to follow the evolution of the secular instability only 
from the linear to weakly nonlinear stages of a monochromatic frequency of gravitational waves in contrast to that 
adopted in this paper. This is because their basic equations are formulated assuming that only one nonaxisymmetric 
mode exists (i.e., assuming that the oscillation frequency of the quadrupole moment is monochromatic). Thus, from 
their simulation, it seems to be difficult to accurately determine the final outcome after the growth of the secular 
instability saturates in which the frequency of quadrupole moment and gravitational waves is not monochromatic 
but gradually changed. On the other hand, their method may be robust to study the early growth of the secular 
instability in which the frequency of gravitational waves is monochromatic. 
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FIG. 3. Snapshots of the density contour curves in the equatorial plane for Ca = 0.3 and e = 1. ?7o « 0.01 in this case. The 
contour curves are drawn for p/pc = O.lj for j = 0.01, 0.1, 1, 2, • • • , 8, 9, 9.5, and 10. Vectors indicate the local velocity field 



1/2 

(v^' , t'-'). and the scale is shown in the upper right-hand corner in units of pc Re- Po denotes the central rotational period of 
the equilibrium configuration given at t = (Po ~ 4.674 in units oi G = Re = Pc = !)■ 
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(a) (b) 

FIG. 5. (a) The density and angular velocity profile along x and y axes at t = (dotted curves) and t = 12.8Po (solid and 
dashed curves) for Ca = 0.3 and e — 1. The solid and dashed curves denote the profiles along x and y axes, respectively, (b) 
The same as (a) but for Ca = 0.334 and t = 24.9Po- 




FIG. 6. Evolution of J and /3 — T/\W\ for Ca ~ 0.3, and e = 1 (solid curve), 0.7 (long dashed curve), and (dashed curve). 
Jo denotes the initial value of J. 
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